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Resolution  and  algorithmic  influences  on  the  baroclinic  pressure  gradient 
in  finite  element-based  hydrodynamic  models 

K.M.  Dresbacka,  C.A.  Blainb,  and  R.L.  Kolara 

aSchool  of  Civil  Engineering  and  Environmental  Science,  University  of  Oklahoma, 

202  W.  Boyd  St.,  Room  334,  Norman,  OK  73019 

b Oceanography  Division,  Code  7322,  Naval  Research  Laboratory, 

Stennis  Space  Center,  MS  39529 

In  3D  shallow  water  models,  the  baroclinic  pressure  gradient  (BPG)  can  become  un¬ 
stable  or  physically  unrealistic  in  areas  of  rapidly  changing  topography  and/or  density 
fields.  By  using  a  2D  finite  element  x-z  model,  whose  formulation  uses  the  generalized 
wave  continuity  equation,  we  assess  the  impact  on  accuracy  of  three  methods  used  to 
compute  the  BPG:  the  sigma  coordinate  method,  the  2— coordinate  (level  coordinate) 
method  and  a  hybrid  method  that  switches  from  sigma  coordinates  to  2— coordinates  at 
a  prescribed  depth.  Resolution  studies  then  look  at  horizontal  and  vertical  resolution 
independently,  and  the  interplay  of  the  horizontal  and  vertical  resolutions  for  these  three 
methods  of  computing  the  BPG.  Numerical  experiments  are  carried  out  on  several  do¬ 
mains,  from  constant  to  rapidly  changing  bathymetry,  with  two  different  density  fields, 
which  vary  only  in  the  horizontal  or  vertical  directions.  Results  thus  far  indicate  that  the 
2— coordinate  method  provides  the  least  amount  of  error  in  the  solution. 

1.  INTRODUCTION 

In  areas  where  the  topography  changes  rapidly,  such  as  a  seamount  or  continental  rise 
region,  many  three-dimensional  hydrodynamic  models  have  problems  computing  a  stable 
and  realistic  BPG.  The  main  topic  of  this  paper  is  the  calculation  of  the  BPG  term.  The 
motivation  for  this  paper  stems  from  anomalous  results  observed  by  Blain  [4]  in  baroclinic 
Arabian  Gulf  simulations  using  a  wave-continuity  based  FE  model.  Figure  1  illustrates 
the  problem,  where  the  source  of  error  was  identified  as  an  unrealistic  (and  unstable) 
BPG  computed  by  the  model  in  a  region  of  steep  bathymetry  and  density  gradients. 

Four  common  coordinate  systems  are  utilized  in  ocean  models:  sigma  coordinates, 
which  are  terrain-following;  2— coordinates  (also  called  level  coordinates),  which  follow 
a  fixed  depth  [3];  isopycnal,  which  follow  lines  of  constant  density;  and  hybrid,  which 
combines  the  sigma  and  z— coordinates.  Advantages  and  disadvantages  exist  for  all  of 
these  coordinate  systems,  which  several  investigators  have  mentioned  in  their  studies 
(e.g.,  [3,6,16]).  To  summarize,  sigma  coordinates  provide  more  resolution  in  the  shallower 
parts  of  the  domain  and  capture  the  bottom  and  free  moving  surfaces,  thus  allowing 
the  boundary  conditions  to  be  implemented  easily.  The  disadvantage  of  sigma  coordinate 
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Figure  1.  Instabilities  caused  by  errors  in  the  calculation  of  the  BPG  during  simulations 
of  the  Arabian  Gulf  [4] . 


system  involves  the  “hydrostatic  inconsistency”  condition,  first  discussed  in  the  context  of 
oceanography  models  by  Haney  [11].  This  condition  indicates  that  in  areas  of  steep  sloping 
topography,  there  needs  to  be  an  accurate  horizontal  resolution  in  order  to  obtain  stable 
and  realistic  BPG  results.  However,  in  some  cases  the  amount  of  horizontal  resolution 
needed  to  produce  accurate  results  lead  to  high  computational  costs.  If  the  “hydrostatic 
inconsistency”  condition  is  not  met,  then  spurious  modes  tend  to  be  introduced  into  the 
solution  through  truncation  errors  obtained  from  the  transformation  of  the  BPG  term  to 
the  sigma  coordinate  system.  Thus,  in  the  areas  of  sloping  bathymetry,  large  truncation 
errors  can  mask  the  true  BPG  [6].  Suggestions  from  some  researchers  have  reduced  these 
errors,  but  the  problem  has  not  been  completely  solved.  In  some  models,  this  problem  has 
been  reduced  by  subtracting  a  mean  vertical  density  gradient  or  an  area-averaged  density 
from  the  initial  density  field  [16-18]. 

As  for  z— coordinates,  they  do  not  suffer  from  problems  with  the  coordinate  transfor¬ 
mation.  However,  the  disadvantage  of  the  z— coordinate  system  is  its  inability  to  properly 
resolve  the  flow  around  the  bottom  topography  in  areas  of  sloping  bathymetry  (“stair¬ 
step”  resolution),  and  the  correct  flow  at  the  surface  is  often  not  captured  [6].  To  obtain 
an  accurate  BPG  where  sloping  bathymetries  come  into  play  for  ^-coordinates,  several 
researchers  suggest  using  extrapolation  techniques,  e.g.  Beckmann  and  Haidvogel  utilized 
a  Chebyshev  polynomial  [2]. 

Another  method  that  has  been  used  in  global  ocean  models  is  isopycnal  coordinates, 
which  follow  lines  of  constant  density  [10].  This  type  of  coordinate  system  does  well  in 
the  deeper  parts  of  the  ocean  because  the  density  profile  tends  to  be  stably  stratified. 
However,  this  coordinate  system  does  not  do  well  in  shallower  parts  of  the  ocean  due  to 
the  mixed  stratified  nature  and  due  to  the  mixing  and  advection  processes  that  tend  to 
be  dominant  in  this  part  of  the  ocean  [6].  It  also  has  the  same  “stair-step”  resolution 
problem  as  the  z— coordinate  method  because  the  lines  of  constant  density  do  not  follow 
the  topography  changes  [5]. 
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Hybrid  methods  have  been  suggested  that  take  advantage  of  the  strengths  of  sigma 
and  2— coordinate,  which  have  been  used  successfully  in  several  models  (e.g,  [6,16]).  The 
degree  of  hybridization  between  the  two  coordinate  systems  and  the  technique  of  the 
hybridization  can  both  vary  in  the  model.  For  example,  Beckers  [1]  examined  a  hybrid 
scheme  that  used  only  one  z— coordinate  (fixed)  with  sigma  coordinates  above  and  below 
it.  Also,  Spall  and  Robinson  [20]  analyzed  a  hybrid  scheme  that  used  2 — coordinates  in 
the  upper  layers  and  sigma  coordinates  in  the  bottom  layers.  Another  hybrid  scheme, 
which  is  used  by  NCOM  (Naval  Coastal  Ocean  Model)  [16],  applies  sigma  coordinates  in 
the  upper  layers  and  the  ^-coordinates  in  the  bottom  layers.  We  note  that  other  types  of 
hybrid  models  have  been  developed,  such  as  HYCOM  [7],  which  switches  from  isopycnal 
in  deep  water  to  sigma  or  2— coordinates  in  coastal  areas. 

Several  researchers  have  investigated  the  unstable  or  unrealistic  results  of  the  BPG 
term  in  the  context  of  finite  difference  models  (e.g.,  [6,11,18]),  however,  only  a  few  studies 
have  been  done  in  the  context  of  finite  element  or  unstructured  grid  models.  Using  a 
finite  element  model,  Walters  and  Foreman  looked  at  the  influence  of  resolution  on  the 
velocity  field  using  sigma  coordinates,  first  varying  the  horizontal  resolution  for  a  fixed 
vertical  resolution  and  then  vice  versa  [21].  They  determined  that  the  sigma  coordinate 
system  produced  either  second-  or  first-order  accurate  solutions,  depending  on  the  density 
profile,  for  the  continental  shelf  region.  From  their  studies,  they  indicated  that  the  sigma 
coordinate  system  studied  should  be  replaced  with  either  z— coordinates  or  the  density 
field  should  be  post-processed  using  a  density  gradient. 

Fortunato  and  Baptista  evaluated  all  of  the  horizontal  gradients  in  the  momentum 
equation  in  either  sigma  coordinates  or  ^-coordinates  in  a  2D  barotropic  and  baroclinic 
(diagnostic)  model  [9] .  They  determined  that  evaluating  all  of  the  horizontal  gradients  in 
the  sigma  coordinate  system  provided  the  best  approach  in  most  cases;  however,  in  certain 
cases  the  z— coordinates  proved  to  be  better,  in  particular  for  the  case  study  presented  in 
Walters  and  Foreman  [21],  They  also  provided  some  steps  to  obtain  the  proper  horizontal 
resolution  for  a  sigma  coordinate  model  near  steep  bathymetry  gradients  [9]. 

Herein,  we  build  on  this  earlier  work  and  investigate  BPG  calculations  using  several  of 
these  coordinate  systems  in  a  finite  element  model.  This  study  will  only  look  at  methods 
to  calculate  the  BPG  term,  while  all  other  horizontal  gradients  in  the  momentum  equation 
will  utilize  sigma  coordinates.  In  this  respect,  it  differs  from  the  work  done  by  Fortunato 
and  Baptista  [9].  The  study  extends  the  work  done  by  Walters  and  Foreman  [21]  by 
looking  at  other  methods  of  calculating  the  BPG  term. 

2.  BACKGROUND  OF  THE  MODEL 

The  model  utilized  in  this  study  is  a  2D  laterally-averaged  shallow  water  model  that 
uses  the  finite  element  method;  it  follows  the  same  development  steps  as  the  3D  ADCIRC 
model  [13].  It  employs  a  mode  splitting  scheme  in  which  the  external  mode  solves  a  ID 
(depth-averaged)  continuity  equation  for  the  elevation  field  and  a  2D  (x-z)  momentum 
equation  to  resolve  the  velocity  field.  The  depth-averaged  velocity  values  utilized  in  the 
continuity  equation  are  obtained  from  the  integration  of  the  2D  momentum  equation 
results.  We  replace  the  continuity  equation  with  the  generalized  wave  continuity  equation 
to  eliminate  the  spurious  modes  that  occur  with  primitive  finite  element  models  (e.g., 
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[12,13,15]).  The  generalized  wave  continuity  equation  is  as  follows: 


*1+Tdc_ 

dt2  °dx  dx 


d_ 

dx 


d(HUU)  dt  „ d\HU ) 

-^  +  gHi-El-^  +  T*D-ToHU-Bx-Dx 


(1) 


=  0, 


where  t  is  surface  elevation  above  a  datum,  H  is  the  total  fluid  depth,  U  is  the  depth- 
averaged  velocity,  g  is  gravity,  Dx  is  the  depth-integrated  momentum  dispersion  (momen¬ 
tum  transfer  due  to  a  non-uniform  velocity  profile),  Bx  is  the  depth-integrated  baroclinic 
forcing,  to  is  a  numerical  parameter  that  allows  either  a  pure  wave  form  of  the  equation 
when  the  parameter  is  small  or  the  primitive  form  of  the  continuity  equation  if  the  pa¬ 
rameter  is  large,  and  r3D  —  Ksiipub,  where  Ksiip  is  the  linear  slip  coefficient  and  ub  is  the 
velocity  at  the  bottom  boundary.  This  code  does  not  include  the  atmospheric  forcing  or 
the  Coriolis  forcing  terms.  We  currently  utilize  a  constant  eddy  viscosity  coefficient,  Ei, 
in  the  horizontal  direction. 

The  model  employs  the  non-conservative  form  of  the  momentum  equation,  which  is  as 
follows: 


du  du  du  d(  d 

M+uTx  +  wTz=-9di  +  a-z 


—  bx  +  Tflx , 


(2) 


where  Tsx/p0  =  Ezdu/dx  is  the  vertical  stress  gradient  with  an  eddy  viscosity  parameter¬ 
ization,  mx  =  d/dx(Eid(u)/dx)  is  the  lateral  stress  gradient  also  with  an  eddy  viscosity 
parameterization  and 


d  [c{p~  Po ) 


Po 


dz 


is  the  BPG.  To  evaluate  Equation  2,  we  mapped  the  terms  onto  a  sigma  coordinate  system, 
where  a  ranges  from  a  at  the  surface  (a  =  1)  to  b  at  the  bottom  (b  =  - 1).  The  equation 
in  the  sigma  coordinate  system  is  as  follows: 


(4) 


where  mx  =  d/dxa(Eidu/dxa)  and  wa  incorporates  terms  from  the  variable  transforma¬ 
tion  (see  [13]  for  more  details).  The  evaluation  of  the  mx  term  occurs  along  the  stretched 
surfaces  directly  with  no  coordinate  transformation. 

These  equations  use  C°  linear  finite  elements  for  the  spatial  discretization  with  the  exact 
quadrature  rules.  For  the  temporal  discretization,  a  three  time-level  scheme  centered  at  k 
is  used  in  Equation  (1)  and  a  two  time-level  scheme  centered  at  k+1/2  is  used  in  Equation 
(4).  Except  for  the  BPG,  all  the  horizontal  derivatives  utilize  the  sigma  coordinate  system; 
the  BPG  term  uses  different  coordinate  systems  as  described  below. 

In  this  work,  we  have  isolated  the  BPG  term  (bx  in  Equation  (4))  and  only  implement 
this  term  differently  in  accordance  with  the  different  coordinate  systems.  Initially,  we 
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calculate  a  buoyancy  term  (Equation  (5)),  which  then  is  used  to  evaluate  the  horizontal 
gradient: 

bduyancy  =  g  f  — — —dz.  /'5'\ 

Jz  po  v  ' 

Then  in  2:— coordinates,  the  BPG  is  given  as 

bx  =  J^~(bouyancy}'  (6) 

while  in  the  sigma  coordinate  system,  the  BPG  is  given  as 
,  _  d  dz  d{bouyancy ) 

x  dx„  dx  dz  ’  (7) 

In  this  study,  we  utilize  both  of  these  equations  along  with  a  combination  of  the  two  to 
develop  the  hybrid  scheme.  The  hybrid  scheme  employed  in  our  study  follows  the  method 
used  in  NCOM  [16],  which  switches  from  sigma  coordinates  to  ^-coordinates  as  the 
depth  increases.  In  the  z— coordinate  BPG  implementation,  we  utilize  the  framework  of 
sigma  coordinates  for  all  terms  except  the  BPG,  which  computes  horizontal  gradients  by 
interpolating  the  values  between  adjacent  vertical  nodes  [2,8].  Near  bottom  boundaries,  a 
linear  extrapolation  technique  is  used  in  regions  where  gradients  based  on  z— coordinates 
“run  into  the  ground” .  Also  note  that  in  this  study,  all  test  cases  were  conducted  utilizing 
the  diagnostic  baroclinic  mode  of  the  model. 

3.  NUMERICAL  EXPERIMENTS 

Three  domains  and  two  different  density  profiles  serve  as  the  test  cases  for  this  study. 
In  order  to  evaluate  each  of  the  methods  for  calculating  the  BPG  term,  we  compare  the 
error  behavior  of  the  BPG,  the  horizontal  velocity  field,  and  the  vertical  velocity  field.  We 
utilize  a  L2  norm  to  determine  the  errors  at  several  stations.  Point  comparisons  at  twelve 
locations  are  made  by  linearly  interpolating  results  from  neighboring  nodes.  These  twelve 
errors  are  then  averaged  to  produce  one  value  for  each  grid  resolution.  In  all  of  the  test 
cases,  we  evaluated  the  horizontal  and  vertical  resolutions  independently,  and  then  the 
interplay  between  both  resolutions.  In  each  test  case,  we  compared  the  results  to  a  “true 
solution”,  which  is  obtained  either  from  an  analytical  solution  or  from  grid  refinement 
studies.  We  indicate  below  what  the  “true  solution”  is  for  each  test  case. 

3.1.  Test  case  1 

In  this  test  case,  a  48-/cm  domain  has  a  constant  bottom  slope  between  10  m  at  the 
shallow  end  to  100  m  at  the  deep  end  (approximately  a  2%  slope),  and  the  density 
profile  varies  linearly  from  1026  kg/m 3  (shallow)  to  1028  kg/m3  (deep)  in  the  horizontal 
direction,  with  no  variation  in  the  vertical  direction.  Boundary  conditions  were  no-flux 
land  boundaries  on  both  sides  of  the  domain.  All  nonlinear  terms  were  evaluated  in  the 
equations.  The  eddy  viscosity  parameters  in  both  the  lateral  and  vertical  directions  are 
kept  constant  at  0  and  0.051  ra2/s,  respectively,  and  bottom  friction  is  evaluated  with  a 
linear  slip  condition  that  utilizes  a  Ksiip  value  of  0.05  m/s.  In  this  test  case,  the  GWC 
equation  numerical  parameter,  r0,  was  set  to  0.001  sec-1  and  a  time  step  of  0.1  sec  was 
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Figure  2.  Bathymetry  and  density  pro¬ 
files  for  Test  Case  2.  Density  values,  which 
vary  from  1000  to  1002  kg/m3,  are  shown 
changing  from  blue  in  lighter  water  to  red 
in  the  heavier  waters.  Bathymetry  ranges 
from  10  to  600  m. 


Figure  3.  Bathymetry  and  density  pro¬ 
files  for  Test  Case  3.  Density  values,  which 
vary  from  1000  to  1002  kg/m3,  are  shown 
changing  from  blue  in  lighter  water  to  red 
in  the  heavier  water.  Bathymetry  ranges 
from  10  to  300  m. 


utilized  for  a  simulation  of  a  day.  Results  were  recorded  90  times  over  the  course  of  the 
simulation. 

A  “true  solution”  for  this  test  case  was  obtained  by  refining  the  grid  in  the  horizontal 
direction  until  the  L 2  error  changes  were  within  machine  accuracy  and  had  reached  con¬ 
vergence,  which  occurred  with  a  constant  nodal  spacing  of  approximately  180  m.  This 
served  as  the  “true  solution”  for  the  horizontal  resolution  study.  Similarly  for  the  vertical 
resolution  study,  we  refined  the  grid  in  the  vertical  until  the  L2  error  changes  showed 
that  the  solution  had  converged,  which  occurred  with  129  nodes  in  the  vertical.  For  the 
interplay  study,  a  nodal  spacing  of  approximately  180  to  in  the  horizontal  with  129  nodes 
in  the  vertical  was  used  for  the  “true  solution.” 

3.2.  Test  case  2 

The  second  test  case  is  taken  from  both  Walters  and  Foreman  [21]  and  Fortunato  and 
Baptista  [9] ,  which  provide  a  test  case  that  mimics  the  shelf  break  region.  The  bathymetry 
varies  linearly  in  three  different  areas  along  a  50  km  slice  as  shown  in  Figure  2.  Density 
varies  only  in  the  vertical  and  depends  on  depth,  as  shown  in  Figure  2.  The  density  field 
above  100  to  comes  from  the  following  relationship:  1001  -  cos(0.0l7rz),  while  below  100 
to,  a  constant  density  value  of  1002  kg/m3  is  used.  Boundary  conditions  were  no-flux 
land  boundaries  on  both  sides  of  the  domain.  All  nonlinear  terms  were  evaluated  in  the 
equations.  The  eddy  viscosity  parameters  in  both  the  lateral  and  vertical  directions  are 
kept  constant  at  0  and  0.051  m2/s,  respectively,  and  bottom  friction  is  evaluated  with  a 
linear  slip  condition  that  utilizes  a  Ksup  value  of  0.001  m/s.  In  this  test  case,  the  GWC 
equation  numerical  parameter,  r0,  was  set  to  0.001  sec-1  and  a  time  step  of  0.1  sec  was 
utilized  for  a  simulation  of  a  day.  Results  were  recorded  90  times  over  the  course  of  the 
simulation. 

For  this  test  case,  we  compared  the  results  to  an  analytical  solution.  As  noted  by 
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Walters  and  Foreman  [21]  and  Fortunato  and  Baptista  [9],  the  analytical  solution  for  this 
test  case  is  zero  because  there  are  no  boundary  forcings  and  the  density  varies  only  in  the 
vertical  direction  (stable  stratification)  and  not  horizontally. 

3.3.  Test  case  3 

In  this  test  case,  we  analyzed  a  different  bathymetry  with  the  same  density  profile 
as  the  second  test  case.  For  the  bathymetry,  we  developed  a  domain  that  includes  a 
seamount,  along  with  a  change  in  topography  that  mimics  the  continental  shelf  region  (see 
Figure  3).  Boundary  conditions  were  a  no-flux  land  boundary  and  elevation  boundary 
of  zero  on  the  opposite  side  of  the  domain.  All  nonlinear  terms  were  evaluated  in  the 
equations.  Eddy  viscosity  parameters  in  both  the  lateral  and  vertical  directions  are  kept 
constant  at  0  and  0.051  m2/s,  respectively.  The  bottom  friction  was  evaluated  with  a 
linear  slip  condition  that  utilized  a  KsHp  value  of  0.001  m/s.  Here  again,  the  GWC 
equation  numerical  parameter,  r0,  was  set  to  0.001  sec-1  and  a  time  step  of  0.1  sec  was 
utilized  for  a  simulation  of  a  day.  Results  were  recorded  90  times  over  the  course  of  the 
simulation. 

As  with  the  second  test  case,  we  compared  the  results  from  this  test  case  to  an  analytical 
solution  of  zero,  because  there  are  no  boundary  forcings  and  the  density  varies  only  in 
the  vertical  direction  (stable  stratification)  and  not  horizontally. 

4.  EXPERIMENTAL  RESULTS 

4.1.  Test  case  1 

The  results  of  the  horizontal  resolution  study  for  the  linearly  sloping  bathymetry  showed 
that  the  only  significant  difference  between  the  three  methods  for  calculating  the  BPG  was 
in  the  value  of  the  BPG  itself.  For  all  the  refinements  in  the  horizontal,  the  z— coordinate 
method  exhibits  a  continual  decrease  in  the  BPG  error;  however,  sigma  coordinates  and 
the  hybrid  scheme  show  that  the  BPG  error  decreases  rapidly  until  the  nodal  spacing  is 
approximately  6000  m  and  then  decrease  slowly  for  the  more  refined  grids.  The  averaged 
Li  errors  for  the  horizontal  and  vertical  velocity  fields  do  not  show  any  appreciable  changes 
based  on  the  different  methods  for  calculating  the  BPG.  However,  when  looking  at  time 
series  of  horizontal  velocity  results,  they  indicate  that  there  are  subtle  changes  between 
the  methods,  yet  these  changes  are  masked  when  computing  the  average  error  over  space 
and  time.  Similar  results  were  found  for  the  norm. 

Results  of  the  vertical  resolution  study  are  nearly  identical  for  the  different  ways  of 
evaluating  the  BPG,  as  was  expected  since  the  variation  of  the  density  field  is  only  based 
in  the  horizontal  direction.  Similarly,  for  this  model  problem,  there  are  no  observable 
changes  between  the  different  methods  when  looking  at  the  interplay  of  the  horizontal 
and  vertical  resolution. 

4.2.  Test  case  2 

Figures  4  and  5  illustrate  the  results  of  the  horizontal  resolution  study  for  the  second 
test  case  for  the  BPG  and  horizontal  velocity.  As  can  be  seen,  the  evaluation  of  the  BPG 
with  coordinates  produced  the  best  results,  while  sigma  coordinates  produced  greater 
errors.  This  greater  sigma  error  is  expected  because  sigma  coordinates  are  more  prone  to 
errors  in  evaluating  the  BPG  when  the  density  field  only  varies  in  the  vertical  direction. 
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Figure  4.  BPG  error  for  Test  Case  2  for 
horizontal  resolution  (vertical  resolution 
held  constant)  for  sigma  coordinates  (short 
dashes),  hybrid  scheme  (long  dashes)  and 
coordinates  (solid  line). 
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Figure  5.  Horizontal  velocity  error  for  Test 
Case  2  for  horizontal  resolution  (vertical 
resolution  held  constant)  for  sigma  coordi¬ 
nates  (short  dashes),  hybrid  scheme  (long 
dashes)  and  2— coordinates  (solid  line). 


The  stretching  of  the  coordinate  system  in  the  vertical  causes  there  to  be  two  different 
density  values  between  the  two  adjacent  sigma  nodes  [9].  The  hybrid  scheme  shows  that 
as  refinement  increases  in  the  horizontal  direction  the  error  decreases  and  approaches 
the  errors  that  are  found  using  ^-coordinates.  The  hybrid  scheme  in  this  case  switches 
from  the  sigma  coordinates  to  z— coordinates  in  the  upper  layer  (at  0  =  —25  m),  where 
there  are  changes  in  the  density.  These  results  follow  that  of  Fortunato  and  Baptista 
[9]  in  indicating  that  ^-coordinates  provide  the  best  solution  to  this  test  case;  however, 
the  hybrid  scheme  shows  promising  results.  Also  note  that  errors  in  the  BPG  produced 
corresponding  errors  in  the  horizontal  velocity  fields  (as  shown  in  Figure  5).  Similar  error 
behavior  is  seen  in  the  vertical  velocity  field  (not  shown). 

We  also  evaluated  the  influence  of  the  vertical  resolution  on  the  BPG  and  velocity  fields. 
Figures  6  and  7  show  the  results  for  the  BPG  and  horizontal  velocity.  As  can  be  seen,  the 
errors  for  coordinates  and  hybrid  scheme  decrease  when  vertical  resolution  is  added; 
however,  we  find  that,  as  further  vertical  resolution  is  added,  the  sigma  coordinate  error 
starts  to  separate  further  from  the  other  two  BPG  calculation  methods  and  appears  to 
reach  an  asymptotic  value.  Similar  error  behavior  is  seen  in  the  vertical  velocity  field  (not 
shown). 

Next,  we  examined  the  interplay  of  the  horizontal  and  vertical  resolution  and  their  effect 
on  the  BPG  and  velocity  errors.  Results  (not  shown)  indicate  that  the  hybrid  scheme  and 
2— coordinates  had  similar  errors  as  refinement  occurred  in  both  horizontal  and  vertical 
directions,  while  the  sigma  coordinate  error  is  higher  by  approximately  one-half  log  cycle 
than  the  other  methods. 

Lastly  for  this  test  case,  we  evaluated  the  placement  of  the  depth  at  which  the  hybrid 
method  switches  coordinate  systems.  Results  show  that  the  depth  should  be  between  20 
m  and  40  m  to  provide  the  lowest  error.  This  depth  range  corresponds  to  the  region  above 
the  rapid  change  in  the  density  field. 


1763 


number  of  nodes 


Figure  6.  BPG  error  for  Test  Case  2  for 
vertical  resolution  (horizontal  resolution 
held  constant)  for  sigma  coordinates  (short 
dashes),  hybrid  scheme  (long  dashes)  and 
^-coordinates  (solid  line). 
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Figure  7.  Horizontal  velocity  error  for  Test 
Case  2  for  vertical  resolution  (horizontal 
resolution  held  constant)  for  sigma  coordi¬ 
nates  (short  dashes),  hybrid  scheme  (long 
dashes)  and  z— coordinates  (solid  line). 


4.3.  Test  case  3 

For  this  test  case,  the  horizontal  resolution  study  showed  similar  results  as  the  second 
test  case  with  the  sigma  coordinate  error  being  higher  than  the  other  two  BPG  calculation 
methods  (Figures  8  and  9).  BPG  results  indicate  that  errors  in  the  hybrid  scheme  (which 
switches  from  sigma  coordinates  to  2: -coordinates  at  2  =  -25  m)  coincide  with  those 
of  the  2:— coordinate  method  as  the  grid  is  refined,  while  the  sigma  coordinate  error  is 
approximately  one-half  log  cycle  higher.  These  BPG  errors  translate  to  similar  errors  in 
the  horizontal  velocity  for  all  the  methods  (cf.  Figures  8  and  9).  We  found  similar  results 
in  the  error  behavior  for  the  vertical  velocity  (not  shown). 

Figures  10  and  11  indicate  how  vertical  resolution  affects  the  error  in  calculating  the 
BPG  and  horizontal  velocity.  The  behavior  of  the  three  BPG  calculation  methods  ex¬ 
hibits  the  same  trends  as  in  the  previous  test  case  with  ^-coordinates  and  the  hybrid 
scheme  having  less  error  than  sigma  coordinates.  Also,  the  sigma  coordinate  error  tends 
toward  an  asymptote  for  both  the  BPG  and  horizontal  velocity.  Results  for  the  vertical 
velocity  show  similar  error  behavior.  In  comparing  the  vertical  resolution  results  from  the 
second  test  case  to  the  results  from  this  test,  one  noticeable  difference  is  that  the  sigma 
coordinate  error  tends  to  separate  itself  more  from  other  two  BPG  calculation  methods. 
More  specifically  in  Test  Case  2,  we  find  that  as  the  vertical  resolution  is  relaxed  the 
errors  approach  the  same  value  for  each  of  the  BPG  calculation  method  but  not  in  Test 
Case  3  (cf.  Figures  7  and  11). 

Next,  we  analyzed  the  interplay  of  the  horizontal  and  vertical  resolution.  BPG  errors 
(not  shown)  indicate  that  in  the  grids  of  lower  resolution  (both  horizontal  and  vertical 
directions)  the  three  methods  of  calculating  the  BPG  approach  similar  values;  however, 
as  more  resolution  is  added  in  both  directions  the  hybrid  scheme  and  2: -coordinates  show 
similar  errors,  while  the  sigma  coordinate  error  is  higher  than  the  other  methods.  In 
the  horizontal  velocity  results,  we  find  that  at  the  lower  resolutions  the  sigma  coordinate 
error  is  slightly  lower  than  the  other  two  methods  of  calculating  the  BPG;  however,  as 
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Figure  8.  BPG  error  for  Test  Case  3  for 
horizontal  resolution  (vertical  resolution 
held  constant)  for  sigma  coordinates  (short 
dashes),  hybrid  scheme  (long  dashes)  and 
2— coordinates  (solid  line). 


Figure  9.  Horizontal  velocity  error  for  Test 
Case  3  for  horizontal  resolution  (vertical 
resolution  held  constant)  for  sigma  coordi¬ 
nates  (short  dashes),  hybrid  scheme  (long 
dashes)  and  2— coordinates  (solid  line). 


Figure  10.  BPG  error  for  Test  Case  3  for 
vertical  resolution  (horizontal  resolution 
held  constant)  for  sigma  coordinates  (short 
dashes),  hybrid  scheme  (long  dashes)  and 
z— coordinates  (solid  line). 
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Figure  11.  Horizontal  velocity  error  for 
Test  Case  3  for  vertical  resolution  (hori¬ 
zontal  resolution  held  constant)  for  sigma 
coordinates  (short  dashes),  hybrid  scheme 
(long  dashes)  and  2— coordinates  (solid 
line). 
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grid  refinement  occurs  the  sigma  coordinate  error  is  higher  than  that  of  the  other  two 
methods  of  calculating  the  BPG.  Similar  behavior  is  found  in  the  vertical  velocity  results. 

Lastly  for  this  test  case,  we  also  evaluated  the  placement  of  the  depth  at  which  the 
hybrid  method  switches  coordinate  systems.  Results  indicate  the  depth  should  be  between 
10  m  to  30  m  for  minimal  error.  As  with  Test  Case  2,  this  depth  range  again  corresponds 
to  the  region  above  the  rapid  change  in  the  density  field. 

5.  CONCLUSIONS  AND  FUTURE  WORK 

Herein,  we  present  some  initial  results  from  an  assessment  of  how  the  coordinate  method 
used  to  calculate  the  BPG  and  resolution  impact  simulation  results.  This  study  uses  a  2D 
laterally-averaged  model  to  investigate  these  changes.  Evidence  thus  far  indicates  that 
the  z— coordinate  method  for  calculating  the  BPG  provides  the  best  results.  However,  the 
three  test  cases  used  in  this  paper  are  more  favorable  toward  the  z— coordinate  system, 
so  in  future  work,  we  will  look  at  other  test  cases  in  which  the  density  and  bathymetry 
profiles  mimic  more  “real”  world  applications,  as  in  the  Arabian  Gulf  problem  presented 
in  the  introduction.  We  will  also  evaluate  variable  grids  with  these  methods  for  calculating 
the  BPG.  Finally,  we  plan  to  look  at  alternative  error  measures,  such  as  the  Cumulative 
Area  Faction  Error  (CAFE)  plots  [14]  or  the  grid  convergence  index  (GCI)  [19],  which 
may  be  less  likely  to  mask  errors  through  averaging  procedures. 

ACKNOWLEDGMENTS 

Financial  support  for  this  research  was  provided,  in  part,  by  a  United  States  Depart¬ 
ment  of  Education  GAANN  Fellowship,  the  Department  of  Defense  under  contract  ONR 
N000 14-02- 1-0651,  and  the  University  of  Oklahoma.  Any  opinions,  findings  and  conclu¬ 
sions  or  recommendations  expressed  in  this  material  are  those  of  the  authors  and  do  not 
necessarily  reflect  those  of  the  funding  agencies. 

REFERENCES 

1.  J.M.  Beckers,  “Applications  of  the  GHER  3D  General  Circulation  Model  to  the  West¬ 
ern  Mediterranean”,  Journal  of  Marine  Systems,  1,  315-322  (1991). 

2.  A.  Beckman,  D.B.  Haidvogel,  “Numerical  Simulation  of  Flow  Around  a  Tall  Isolated 
Seamount.  Part  1:  Problem  Formulation  and  Model  Accuracy”,  Journal  of  Physical 
Oceanography,  23,  1736-1753  (1993). 

3.  M.D.J.P.  Bij velds,  J.A.Th.M.  van  Kester,  G.S.  Stelling,  “A  Comparison  of  Two  3D 
Shallow  Water  Models  Using  Sigma-coordinates  and  z— coordinates  in  the  Vertical 
Direction” ,  Estaurine  and  Coastal  Modeling,  Proceedings  of  the  Sixth  International 
Conference,  M.L.  Spaulding  and  H.L.  Butler  (eds.),  American  Society  of  Civil  Engi¬ 
neering,  130-147  (2000). 

4.  C.A.  Blain,  personal  communication  2001. 

5.  R.  Bleck,  “Finite  Difference  Equations  in  Generalized  Vertical  Coordinates.  Part 
1:  Total  Energy  Conservation”,  Contributions  in  Atmospheric  Physics ,  51,  360-372 
(1978). 

6.  H.  Burchard,  O.  Petersen,  “Hybridization  Between  s-  and  z— Co-ordinates  for  Improv- 


1766 


ing  the  Internal  Pressure  Gradient  Calculation  in  Marine  Models  with  Steep  Bottom 
Slopes”,  International  Journal  for  Numerical  Methods  in  Fluids,  25, 1003-1023  (1997). 

7.  E.P.  Chassignet,  L.T.  Smith,  G.R.  Halliwell,  R.  Bleck,  “North  AtlanticSimulations 
with  the  HYbrid  Coordinate  Ocean  Model  (HYCOM):  Impact  of  the  Vertical  Coordi¬ 
nate  Choice,  Reference  Pressure,  and  Thermobaricity” ,  Journal  of  Physical  Oceanog¬ 
raphy,  33,  2504-2526  2003. 

8.  A.B.  Fortunato,  A.M.  Baptista,  “Localized  Sigma  Coordinates  for  the  Vertical  Struc¬ 
ture  of  Hydrodynamic  Models”,  in  Estaurine  and  Coastal  Modeling,  Proceedings  of 
the  Third  International  Conference,  M.L.  Spaulding  et  al.(eds.),  American  Society  of 
Civil  Engineering,  323-335  (1994). 

9.  A.B.  Fortunato,  A.M  Baptista,  “Evaluation  of  the  Horizontal  Gradients  in  Sigma 
Coordinates  in  Shallow  Water  Models”,  Atmosphere-Ocean,  34(3),  489-514  (1996). 

10.  D.  Haidvogel,  Numerical  Ocean  Circulation  Modeling,  J.N.B.  Bell  (ed.),  Beckman 
Imperial  College  Press,  London,  320  (1999). 

11.  R.L.  Haney,  “On  the  Pressure  Gradient  Force  over  Steep  Topography  in  Sigma  Coor¬ 
dinate  Ocean  Model”,  Journal  of  Physical  Oceanography,  21,  610-619  (1991). 

12.  I.P.E.  Kinnmark,  “The  shallow  water  wave  equations:  formulations,  analysis  and 
application”,  in  Lecture  Notes  in  Engineering,  C.A.  Brebbia,  S.A.  Orsszag  (eds), 
Springer- Verlag,  Berlin,  15,  187  (1986). 

13.  R.A.  Luettich,  J.J.  Westerink,  “Formulation  and  Numerical  Implementation  of  the 
2D/3D  ADCIRC  Finite  Element  Model  Version  43.02”,  internal  report,  54,  (2003). 

14.  R.A.  Luettich,  J.J.  Westerink,  “Continental  Shelf  Scale  Convergence  Studies  with  a 
Tidal  Model”,  in  Quantitative  Skill  Assessment  for  Coastal  Ocean  Models,  47,  D.R. 
Lynch  and  A.M.  Davies(eds.),  American  Geophysical  Union:  Washington,  DC,  349- 
371  (1995). 

15.  D.R.  Lynch,  W.G.  Gray,  “A  Wave  Equation  Model  for  Finite  Element  Tidal  Compu¬ 
tations”,  Computers  and  Fluids,  7(3),  207-228  (1979). 

16.  P.J.  Martin,  “Description  of  the  Navy  Coastal  Ocean  Model  Version  1.0”,  NRL  Report 
No.  NRL/FR/7322-00-9962,  45  (2000). 

17.  G.L.  Mellor,  T.  Ezer,  L.Y.  Oey,  “The  Pressure  Gradient  Conundrum  of  Sigma  Coordi¬ 
nate  Ocean  Models”,  Journal  of  Atmospheric  and  Oceanic  Technology,  11,  1126-1134 
(1994). 

18.  G.L.  Mellor,  L.Y.  Oey,  T.  Ezer,  “Sigma  Coordinate  Pressure  Gradient  Errors  and  the 
Seamount  Problem”,  Journal  of  Atmospheric  and  Oceanic  Technology,  15,  1122-1131 
(1998). 

19.  P.J.  Roache,  Verification  and  Validation  in  Computational  Science  and  Engineering, 
Hermosa  Publishers,  Albuquerque,  New  Mexico,  446  1998. 

20.  M.A.  Spall  and  A.R.  Robinson,  “Regional  Primitive  Equation  Studies  of  the  Gulf 
Stream  Meander  and  Ring  Formation  Region” ,  Journal  of  Physical  Oceanography, 
20,985-1016  (1990). 

21.  R.A.  Walters,  M.G.G.  Foreman,  “A  3D,  Finite  Element  Model  for  Baroclinic  Circu¬ 
lation  on  the  Vancouver  Island  Continental  Shelf”,  Journal  of  Marine  Systems,  3, 
507-518  (1992). 


